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A  COMPUTER  CODE  FOR  HIGH  ALTITUDE  EMP 

I .  Introduction 

The  effects  of  a  nuclear  environment  on  aerospace 
systems  is  an  important  factor  in  systems  analysis.  During 
the  past  few  years  several  students  have  worked  with 
Professor  Bridgman  at  the  Air  Force  Institute  of  Technology 
(AFIT)  on  a  computer  code  to  determine  survivability  of  a 
system  with  known  nuclear  vulnerabilities  from  a  variable 
nuclear  threat.  <!he  AFIT  survivability  code  capabilities 
include  blast,  thermal,  x-ray,  gamma  ray,  and  neutron 
effects.  The  high  altitude  EMP  code  presented  in  this 
report  ii  intended  to  be  used  in  conjuction  with  the  AFIT 
survivability  code. 

The  EMP  (electromagnetic  pulse)  from  a  nuclear  weapon 
is  usually  considered  to  be  a  radiating  electromagnetic 
wave  of  short  duration  containing  many  frequencies. 

However,  the  nuclear  generated  EMP  was  not  studied  seriously 
until  a  considerable  time  after  the  first  nuclear  explosion. 

At  present  there  is  a  significant  amount  of  work  being  done 
to  model  EMP  generation  and  effects.  For  example  the  Air 
Force  Weapons  Laboratory  (AFWL)  and  several  civilian 
companies  under  contract  to  the  USAF  are  working  in  the  field. 

There  are  several  different  types  of  EMP  with  distinc¬ 
tions  made  between  the  mechanisms  which  generate  them. 

Kinsley  (Ref  1)  presents  a  comprehensive  discussion  of 
the  various  types  of  EMP.  For  example,  a  nuclear  burst 
on  the  ground  produces  an  EMP  with  different  characteristics 
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than  those  from  a  high  altitude  burst.  Also,  nuclear  burst 
products  interacting  directly  with  a  system  can  produce  an 
EMP  within  the  system  or  even  within  the  circuits  of  the 
system.  This  report  considers  only  the  EMP  generated  by 
high  altitude  burst  gamma  rays  interacting  with  the  at¬ 
mosphere  . 

The  high  altitude  EMP  code  developed  in  this  report 
is  based  on  the  theory  of  Karzas  and  Latter  (Ref  2). 
Briefly,  the  theory  develops  a  model  for  the  interaction 
of  Cpmpton  electrons  with  the  geomagnetic  field.  The 
Compton  electrons  are  produced  by  prompt  gamma  radiation 
from  the  burst  in  a  reasonably  well  defined  region  in  the 
atmosphere.  Several  simplifications  are  made  before 
arriving  at  the  final  equations. 

Since  several  of  the  simplifications  and  assumptions 
used  are  implicit  in  the  presentation  of  the  theory,  it 
is  appropriate  to  list  them  here.  Only  one  group  of 
monoenergetic  unscattered  gamma  rays  are  considered  to 
produce  Compton  electrons.  Each  gamma  which  interacts 
is  assumed  to  produce  one  and  only  one  Compton  electron 
initially  traveling  precisely  in  the  radial  direction. 

No  angular  distribution  of  Compton  electrons  is  allowed. 

All  Compton  electrons  are  assumed  to  have  the  same  energy. 
Curvature  of  the  Earth's  magnetic  field  is  ignored.  The 
electromagnetic  fields  are  not  self-consistent,  that  is, 
only  the  geomagnetic  field  is  considered  to  affect  the 
motion  of  the  Compton  electrons.  Cascading  of  secondary 
electrons  and  recombination  of  ions  is  ignored.  The  low 
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frequency  portion  of  the  pulse  is  not  considered.  The  Earth 
is  assumed  to  be  flat  and  the  finite  conductivity  of  ground 
is  not  considered.  The  burst  is  assumed  to  be  far  from  the 
absorption  region.  Only  gamma  ray  effects  are  considered. 

Although  the  final  model  is  somewhat  restricted  by 
these  assumptions  and  simplifications,  the  end  result  is 
a  relatively  inexpensive  computer  code  which  gives  a 
peak  value  of  the  electric  field  at  any  target  point  on  or 
above  the  ground,  which  is  an  upper  bound  on  the  actual 
peak  value. 

Section  II  of  this  report  develops  the  theory  and 
derives  the  equations  used  in  the  code.  Section  III 
describes  the  calcul ational  procedures  used  in  the  code. 
Section  IV  presents  a  sample  of  typical  results  and  a 
study  of  input  parameter  variation.  Section  V  is  a 
discussion  of  the  code's  limitations  and  uses,  with 
recommendations  for  possible  improvements.  Appendix  A 
is  a  code  user's  guide.  Appendix  B  is  the  detailed  flow 
charts  for  the  entire  code.  And  finally  Appendix  C  is  a 
listing  of  the  complete  code. 
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II./..  g^eory 

Overview 

The  EMP  problem  is  a  problem  in  classic  electromagnetic 
theory.  A  solution  of  Maxwell's  equations  is  a  solution 
of  the  problem.  In  this  case  it  is  necessary  to  model  the 
current  and  charge  densities  generated  by  the  gamma  rays  in 
the  absorption  region  to  obtain  the  sources  and  conductivity 
needed  to  solve  Maxwell's  equations. 

Expressions  for  the  current  sources  and  conductivity 
are  obtained  in  four  steps.  The  transport  of  the  gammas 
from  the  burst  to  the  absorption  region  is  used  to  obtain 
the  number  density  of  reacting  gammas.  This  result  is  used 
with  the  models  for  the  current  and  charge  densities  to 
obtain  preliminary  expressions.  Then  after  considering 
the  relativistic  notion  of  the  Compton  electrons,  the 
preliminary  expressions  are  transformed  to  spherical  co¬ 
ordinates  . 

After  presenting  Maxwell's  equations  in  a  convenient 
fora,  they  are  transformed  to  spherical  coordinates  and 
retarded  time.  A  high  frequency  approximation  is  then  made 
to  arrive  at  the  final  equations. 

Electron  Current  and  Density 

Ganna  Transport .  Consider  the  geometry  shown  in  Fig.  1. 
The  nuclear  burst  occurs  at  the  origin  at  time,  t  »  0.  The 
gamma  rays  move  to  point  r'  in  tine  t'  and  at  that  point  and 
time  interact  to  create  Compton  electrons.  It  is  assumed 
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that  each  ganna  creates  one  and  only  one  Conpton  electron 
traveling  in  the  radial  direction  with  the  naxinun  Conpton 
recoil  energy. 

The  ganna  ray  emission  rate  can  be  taken  as 

■  F  C») 

where 

N (t)  ■  nunber  of  ganna  rays  enitted 

Y  ■  ganna  ray  yield  of  burst 

E  ■  nean  energy  of  the  ganna  rays 

f(t)  »  tine  dependence  of  the  yield 

and 

/"  f(t)dt  ■  1  (2) 


The  nunber  density  of  gannas,  g(r),  which  interact 
at  a  point,  r,  can  be  taken  as 

i  rT  dr' 

l(r)  -  l 

4wr2X(r) 

where 

X  «  nean  free  path  for  production  of  Conpton 
electrons. 

Electron  Currents  and  Densities .  The  rate  of  product! 
of  primary  (Conpton)  electron  density,  "pri* 

“If1*-  -  g(r)f(t  -  r/c)  (4) 
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Following  the  Karzas-Latter  approach  (Ref  2)  it  is 
assuaed  that  the  electrons  maintain  their  initial  speed, 

Vq,  throughout  their  range,  R,  and  then  abruptly  stop. 

Also,  it  is  assuaed  that  the  secondary  electrons  are  aade 
at  a  unifora  rate  during  the  lifetime,  R/VQ,  of  the  Conpton 
electrons.  Therefore,  the  rate  of  production  of  secondary 
electron  density,  n, . . ,  is 

SOv 


J|I5.C 


■I 


Epri/33ev 

r/v;- 


n 


pri 


«V0 

T 


n 


pri 


(S) 


where 

Epri  ■  the  initial  energy  of  the  Compton  electrons 
R  ■  the  range  of  the  Compton  electrons  in  air 

«  ■  Epri/J5,v 

33ev  ■  average  ionization  energy  per  molecule  for 
air 

Vq  »  the  speed  of  the  Compton  electrons 
R/Vq  ■  the  lifetime  of  the  Compton  electrons 

Now  consider  the  current  resulting  from  the  Conpton 
electrons.  The  differential  current  is  the  charge  times 
the  velocity  times  the  differential  density  of  electrons. 
Hence 


dJc  -  -eV(t-t 'JgCr'JfCt*  -  r'/c)dt' 


(6) 
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where 

♦ 

V(t-t')  ■  velocity  of  the  Conpton  electrons  at 
tine  t  which  were  crested  at  tine  t'. 


Putting  (6)  into  integral  fora  gives 

JC  "  "e  ^t-R/VQ  *(r’)f£t'  -  r'/c)$(t-f)df  (7) 

Now  let 

t»  ■  t  -  t'  (8a) 

t  ■  t  -  (r/c)  (8b) 

X(T')  -  X(t)  •  r  -  r*  (8c) 

Also  note  that 

(r-r')  «  r  or  r*  (9) 

for  distant  explosions  (see  Fig.  1).  So, 

g(r)  -  l(r ' )  (10) 

Using  Eqs  (8) ,  (9) ,  and  (10)  in  Eq  (7)  gives 

Jc  ■  -eg(r)  /q/V°  V(T*)f^T  -  t»  ♦  MilL^dt'  (11) 

Using  similar  arguments, 

npri  "  8(r)  f0  °  f  (T  •  T'  ♦  dT'  C*2) 
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And  putting  Eq  (12)  into  Eq  (5)  yields 
qVQ  X 

".oc  ■  T-  "pri‘T,)dt' 

■  g(r)  ^*/V°  *  (»'  -  T"  ♦  j  dx"l  dt’  (IS) 

Relativistic  Electron  Motion.  Equations  (11),  (12), 
and  (13)  contain  r(T)  and  X(t)  which  are  not  yet  defined  in 
an  easily  obtained  form.  The  equation  of  aotion  for  a 
Compton  electron  is 

Jt  (ayV)  -  -eVXB0  (14) 

where 

a  •  electron  rest  aass 

y  •  [l-(V/c)Vl/2 

Bq  ■  earth's  aagnetic  field  » 

Again  it  is  assuaed  that  VQ  is  constant  throughout  the 
electron's  lifetiae. 

With  at  *  eBQ/ay  Eq  (14)  becoaes 

^  V(t)  -  -V(T)XUxw  (15) 

Breaking  Eq  (15)  into  its  rectangular  coaponents 

dVx 

JT  -  -«Vy  (16a) 
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dV 

3r  ■  “vx 

(16b) 

A  solution 

dV 

jr-' 0 

for  this  set  of  equations  is 

(16c) 

• 

V  ■  Vi  cos  at 
x  ■*' 

(17a) 

V  ■V,  sin  wt 

y  ± 

(17b) 

v,  ■  v« 

(17c) 

where  is  the  initial  velocity  coaponent  perpendicular  to 

♦  ♦ 
Bq  and  V||  is  the  initial  velocity  coaponent  parallel  to  BQ 

and  both  are  constants  with  respect  to  t. 

Transfornation  to  Spherical  Coordinates.  It  is  conveni¬ 
ent  to  transfora  the  above  solution  to  a  spherical  coordinate 
systea  with  its  origin  at  the  burst  point.  The  transforaa~ 
tion  froa  rectangular  to  spherical  coordinates  is 

V  »  V  sin  0  cos  0  ♦  Vw  sin  0  sin  0  ♦  V,  cos  0  (18a) 

r  x  y  * 

Vg  ■  cos  0  cos  0  ♦  Vy  cos  0  sin  ♦  *  vx  cos  0  (18b) 
V.  ■  -V  sin  0  *  V  cos  0  (18c) 

?  *  y 

Without  loss  of  generality  -the  coordinates  can  be  chosen 
such  that  V  lies  in  the  X-Y  plane,  hence  0*0,  and  the 
transfornation  becomes 

( 
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V  ■  V  sin  0  ♦  V.  cos  0 
r  x  * 

(19a) 

V  ■  V„  cos  0  -  V  sin  0 

9  *  z 

i  19b) 

■  V  i 

\ 

(19c) 

Note  that 

\ 

Yl  ■  v0  •  \ 

(20a) 

*11  *  V0  CM  '  \ 

(20b) 

Putting 

Eqs  (17)  and  (20)  into  Eq  (19)  gives 

V  ■  Vn[sin2  0  cos  wr  ♦  cos2  0] 

T  U 

(21a) 

V.[cos  0  sin  0  cos  wt  -  sin  0  cos  0] 

(21b) 

■  VQ[sin  0  sin  wt]  (21c) 

Now  X(t)  can  be  written  as 

X(T)  "  f0  VrdT  *  V0[,in2  6  'Sl£  ♦  T  co,:  0J  <22) 
Equations  (21)  and  (22)  substituted  into  Eq  (11)  give 

R/V 

■  -eg(r)VQ  /  0  f (T)  [cos2  0  ♦  sin2  0  cos  wT’JdT'  (23) 

R/V 

J*  ■  -«g(r)V.  I  0  f  (T)  [sin  0  cos  0  (cos  WT'-l)]dT»  (24) 

O  U  0 

.  r/va 

■  -«g(*,)V0  /Q  0  f  (T)  [sin  0  sin  WT*]dT*  (25) 


11 


GNE/PH/74-1 


( 
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i 


■ 
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where 

T  -  t  -  (1  -  0  cos2  6)T'  ♦  6  sin2  6  (26a) 

with 

6  -  VQ/c  (26b) 

Equations  (23) ,  (24)  ,  (25) ,  and  (26)  provide  the 
Compton  currents  within  the  absorption  region  in  a  fora  which 
can  be  used  in  the  final  field  equations.  In  addition  to 
the  Compton  currents,  an  expression  for  the  conductivity 
within  the  absorption  region  is  needed. 

Equations  (21)  and  (22)  substituted  into  Eq  (13)  give 

1v0  T  »/»„ 

"».C(T)  ‘  T  *(r)  l/0  (27) 

wher. 

T*  -  T'  -  (1  -  0  cos2  0)t«  ♦  8  Sin2  0  5i±J*X"  (28) 

S 

Consider  the  equation  of  motion  for  secondary  electrons. 

*  -*■ 

Neglecting  the  VXBQ  term,  which  i;>  small  compared  to  the 
other  terns  (Ref  2)  it  is 

■  Jjf  ■  -•£  »Vvc  (29) 

where 

v  ■  electron  collision  frequency, 
c 
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These  secondary  electrons  have  velocities  in  the  theraal 
region  and  are  assumed  to  reach  their  maximum  velocity 
instantly.  In  this  case,  Eq  (29)  becomes 

V  -  jj£-  E  (30) 

c 

The  current  source  from  the  secondary  electrons  is 

Jsec  -  -«V(T)n  (x)  -  e2E  (31) 

sec  "c 

or  in  terms  of  conductivity 

Jsec  -  o(t)E  (32) 

where 

niec(t)  2 

«(T)  -  -555 -  «2  (SJ) 

c 

Equations  (32)  and  (33)  provide  the  needed  expressions 
for  the  conductivity. 

Electromagnetic  Fields  from  High  Altitude  Currents 

Maxwell *s  Equations .  Now  that  the  Compton  currents  and 
the  conductivity  due  to  secondary  electrons  have  been  ob - 
tained,  consider  the  field  equations. 

Maxwell's  equations  are 

$xE  -  -  |i  (34a) 
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1  dl 


VxB  ■  M0J  ♦  -y  7t 
c 


(34b) 


7 • E  -  -1 
£0 


$*B  ■  0 


(34c) 


(34d) 


where 


1  ■  total  current  density 


q^  ■  total  charge  density 


Continuity  of  charge  requires 


♦  ?.j  .  o 


(JS) 


It  is  convenient  to  conbine  the  above  equations  into 
equations  containing  only  E  in  one  and  only  f  in  the  other. 


Doing  so  gives 


(’’  ■  7  S?)  !  -•  ^  ^  !V 


(37) 


Transforastion  to  Spherical  Coordinates  and  Retarded 
Tine.  Equations  (36)  and  (37)  will  now  be  transforned  to 
spherical  coordinates  and  retarded  tine.  Consider  the  trans- 
fornation 
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T  ■  t  -  r/c 


(38*) 


r»  ■  r 


(38b) 


0»  •  0 


(38c) 


♦  '  -  ♦ 


(38d) 


This  is  *  spherical  coordinate  systea  where  tiae  is  aeasured 
at  each  radial  point  in  terns  of  the  arrival  of  the  boab 
gaaaa  rays  at  that  point. 

Using  Eq  (38)  it  is  easily  shown  that 


Thus  the  operator 


transforas  to 


a  .  a 


(39) 


a  a  i  a 

SF  *  “  c  ST 


(40) 


a  a 

ST  •  W 


(41) 


a  a 
5?  •  Wr 


(42) 
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and  the  operator 


transforms  to 


♦  A  i  a 
7  -  U  I  |- 
r  c  er 


Similarly,  the  operator 


1  ii 

c2  3t 


transforms  to 


o2  2  3  1  3 

7  '  c  Sr  r  Jr  r 


Equation  (36)  now  becoaas 


r„2  2  1  3  3  I*  3J  l  *  *  1  8,y 

’  c  F  St  J7  rJ  E  ’  "o  St  ‘  *  ur  c  Sr-  t4s) 


and  Eq  (3S)  becomes 


.  .?.J  ♦  G  1 1- 

3t  r  c  3t 


•V’J  *  isr1 


(44) 


Using  Eq  (44)  in  Eq  (43)  and  rearranging  gives 


2+  A  1  4  4  i  ■+■ 

-V  E  ♦  U  -4-  V  •  J  ♦  4-  Vq 
r  ce„ 


*  I? 


[H 


(rE)  ♦  p. 


,(J-OrJr)  ] 


(45) 
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Similarly,  Eq  (37)  becomes 


■»*»*  •&*ir[hh  <'*)] 

'!f[r'VfW]  -° 


(46) 


High  Frequency  Approximation .  Again,  following  the 
Karzas-Latter  model,  note  that  the  variation  of  currents 
with  distance  is  slow  compared  to  variations  with  time  and 
that  the  fields  resulting  from  the  transverse  currents  are 
rapidly  varying  in  character,  as  are  the  currents  themselves. 
Therefore,  only  the  3/3t  terms  are  kept  in  the  transverse  field 
equations.  Since  the  radial  components  do  not  propagate  out¬ 
side  of  the  absorption  region,  they  are  not  considered  further 
in  this  report. 

The  transverse  equations  become 


§7  [I  F  If  (rE0)  *  Ve  ]  *  0 


(47) 


lx  [I  r  I?  (rE*>  +  V*  ]  "  0 


(48) 


I? [I  f  If  (rV  "  r-  J$] " 0 


(49) 


?r[!F?F(rV  ♦!rJe]-# 


(50) 


17 


( 


These  equations  are  called  the  Karzas-Latter  high 
frequency  approximation  for  the  EHP  fields,  and  they  are 


useful  in  the  range  0  <  t  <  100  shakes. 

Integration  with  respect  to  tiae  yields 

I  F  IF  4  V#  *  0  (5l> 

HI?  <*v  4  v*  • 0 

HI?  -  r  4  0  <S5> 

HI?  4  r  •»»  * 0  <S4> 

Recall  that  the  total  current  density  is  <■ 

J  -  Jpri  ♦  JS0C  •  Jc  ♦  o(t)E  (SS) 

so  Eqs  (SI)  and  (52)  become 

H  I?  <rV  4  "oJS  4  ’  0  (S6) 

hi?  <rv 4  v; 4  %o(t)e*  ■ #  (s7> 

With  the  aid  of  a  computer,  it  is  now  possible  to  obtain 
numerical  solutions  for  the  above  equations  which  will  yield 
a  slightly  high  estimate  of  the  peak  value  of  the  EMP  pulse 
resulting  from  a  high  altitude  burst. 
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Below  the  absorption  region  the  Compton  currents  and 
the  conductivity  are  zero.  In  this  case,  Eqs  (56)  and  (57) 
have  the  following  solutions: 

E#  ■  Cj/r  (SB) 

Eg  -  C2/r  (59) 

where  Cj  and  C2  are  deterained  by  the  values  of  Eg,  Eg,  and 
r  at  the  bottom  of  the  absorption  region. 
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III.  Code  Description 

General  Approach 

Equations  (56),  (57),  (58),  and  (59)  were  chosen  as 
the  simplest  ones  to  solve  nunerically.  Of  course,  Eqs 
(24),  (25),  (27),  and  (33)  are  used  to  obtain  the  Coapton 
currents  and  conductivity  needed  to  solve  Eqs  (56)  and  (57). 
The  B  -  field  equations  are  not  solved  since 

E  -  cB  (60) 

can  be  used  to  obtain  B  once  E  is  found.  This  relationship 
is  based  on  the  assuaption  that  the  EMP  pulse  is  a  spherical 
wave  propagating  in  free  space,  below  the  absorption  region. 

The  function  used  for  the  tine  dependence  of  the  weapon 
yield  is  the  one  recoaaended  by  Poaranning  (Ref  3). 

(*♦0)  e*P  (t-t  ) 

f(T)  ■  (1/N)  g  +  a  #Xp  [  (04-B)  (t-Tq)J 
where  N  is  chosen  such  that 

/"f(T)dt  -  1  (62) 

and  a  >  0. 

This  function  rises  like  eaT  for  snail  t,  falls  like  e’^T 
for  large  t,  and  has  a  single  aaxiaua  at  tq. 

Figure  2  presents  a  flow  chart  which  is  descriptive 
of  the  approach  taken  solving  the  equations.  The  top  of 
the  absorption  region  is  assuaed  to  be  at  50  kn  altitude 
and  the  bottoa  of  the  absorption  region  is  assuaed  to  be 


GNE/PH/74-1 


at  20  km  altitude.  Calculations  by  Latter  and  LeLevier 
(  (Ref  4)  indicate  that  20  km  to  SO  km  is  the  altitude  where 

most  of  the  prompt  gamma  ray  energy  is  deposited. 

Figure  3  depicts  the  target  geometry.  The  value  for 
RMIN  is  determined  by  the  intersection  of  the  line  of  sight 
with  the  50  km  altitude.  The  value  for  RMAX  is  determined 
by  the  intersection  of  the  line  of  sight  with  the  20  km 
altitude.  If  the  target  is  in  the  absorption  region  the 
target  altitude  determines  RMAX  for  the  direct  wave  calcula> 
tion.  These  two  values  of  r  are  the  limits  on  the  mesh  in 
the  r  direction.  The  line  of  sight  is  divided  into  the 
desired  number  of  steps  along  r  for  the  integration  on  r 
in  the  absorption  region. 

The  retarded  time  direction  of  the  mesh  is  divided 

( 

into  0.1  shake  steps  up  to  10  shakes  and  then  1.0  shake 
steps  on  up  to  100  shakes.  Calculation  can  be  stopped  at 
any  desired  TMAX  from  10  to  100  shakes*  which  is  the  upper 
limit  of  the  usefulness  of  the  high  frequency  approximation. 

If  the  ground  reflected  wave  is  to  be  calculated* 
the  mirror  image  of  the  target,  below  ground,  is  used 
to  find  the  line  of  sight  from  the  burst  to  the  target. 
(Refer  to  Fig.  3.) 

At  r  ■  RMIN  all  of  the  fields  are  assumed  to  be  zero. 
For  each  x,  equations  (57)  and  (58)  are  integrated  over  r 
from  RMIN  to  RMAX  and  the  value  of  E  at  the  bottom  of  the 
absorption  region  is  stored.  At  each  step  in  r,  equations 
\  (24),  (25),  and  (27)  are  numerically  integrated.  Then 
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equations  (S8)  and  (S9)  are  combined  into 

„  (UMAX) (Ermax)  (63) 

rt.rg« 

to  find  E  at  the  target. 

The  E  array  is  then  searched  to  find  the  peak  value 
before  printing  out  the  results. 

Inputs 

The  code  uses  a  right  handed  Cartesian  coordinate 
system  with  the  ground  in  the  X-Y  plane,  the  vector  in 

A 

the  Y-Z  plane,  and  j  pointing  towards  the  equator.  For 

A 

example,  in  the  northern  hemisphere,  i  is  magnetic  west, 

A  A 

j  is  magnetic  south,  and  k  is  altitude.  The  origin  of 
the  coordinate  system  is  always  at  ground  zero,  directly 
below  the  burst.  Note  that  this  coordinate  system  is  not 
the  same  as  the  Cartesian  systems  used  earlier. 

Refering  to  the  above  coordinate  system  the  target 
coordinates,  (X,Y,Z),  are  read  in  using  units  of  meters. 
If  the  reflected  wave  is  to  be  calculated  the  altitude  is 
read  in  as  a  negative  number,  (X,Y,-Z). 

The  height  of  the  burst  is  read  in  using  units  of 
kilometers.  The  gamma  yield  of  the  burst  is  read  in 
using  units  of  kilotons. 

The  magnitude  of  the  Earth's  magnetic  field  is  read 
in  using  units  of  webers  per  square  meter.  The  dip  angle 
($  in  Fig.  1)  of  the  magnetic  field  is  read  in  using 
units  of  degrees. 
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NDELR,  the  desired  number  of  steps  to  be  used  in  the 
integration  over  r  in  the  absorption  region,  is  read  in  as 
any  integer  in  the  closed  interval  [50,  500]. 

TMAX,  the  retarded  time  where  calculations  are  to  be 
stopped,  is  read  in,  using  units  of  shukes,  as  any  integer 
in  the  closed  interval  [10,  100]. 

Preliminary  Calculations 

Before  starting  the  numerical  integrations,  the  code 
performs  several  preliminary  calculations.  The  input  data 
is  converted  to  MKS  units.  The  reflected  wave  is  used 
whenever  Z  is  greater  than  49  km  or  less  than  0.  The 
target  coordinates  are  transformed  to  a  spherical  coordinate 
system  with  the  burst  at  the  origin  and  the  polar  axis 
parallel  to  The  line  of  sight  intersections  with  the 

absorption  region  are  determined.  And  finally,  the  constant 
angles  required  by  the  code,  0  and  A,  (see  Fig.  1)  are 
calculated. 

Calculation  of  Compton  Currents  and  Conductivity 

The  two  Compton  currents,  Jq  and  are  calculated 
at  each  r,  t  mesh  point  by  numerically  integrating  equations 
(24)  and  (25).  The  step  size  used  is  0.1  times  the  Compton 
lifetime,  R/Vq.  The  integration  itself  is  done  using  the 
4th  order  Runge-Kutta  method  (Ref  5).  It  should  be  noted 
that  both  the  mean  free  path  for  Compton  interaction  and 
the  Compton  lifetime  are  exponentially  scaled  from  sea 
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level  values  using  a  7  km  scale  height.  However,  the 
Compton  lifetime  is  not  allowed  to  be  greater  than  100 
shakes,  since  this  is  the  maximum  time  of  interest. 

Monoenergetic  gammas  of  energy  1.5  Mev  are  assumed. 
The  most  energetic  Compton  electrons  resulting  from  1.5 

O 

Mev  gammas  have  a  speed  of  2.88  (10)  m/sec.  Therefore 
V0  -  2.88  (10) 8  m/sec. 

Since  the  integration  on  t"  in  equation  (27)  is  also 
over  the  Compton  lifetime,  this  integration  is  carried 
out  simultaneously  with  the  Compton  current  integrations. 
Again,  the  4**1  order  Runge-Kutta  method  is  used.  It 
is  broken  into  two  parts,  one  for  -•  <  t'  <  0  and  the 
other  for  0  <  t'  <  t.  In  this  case,  -•  is  defined  to  be 
the  time  when  the  first  gamma  ray  reached  the  top  of  the 
absorption  region,  since  no  secondaries  can  be  produced 
before  that  time. 

The  integration  on  t'  in  equations  (27)  is  also 
broken  into  two  parts,  one  for  -•  <  t'  <  0  and  the  other 
for  0  <  t'  <  t.  In  the  first  case,  integration  is  started 
at  x*  ■  0  and  proceeds  to  x*  ■  -(r-RMIN)/Vg  in  steps  of 
Ax'  »  -Ar/VQ.  In  the  second  case,  integration  is  started 
at  x'  *  0  and  proceeds  to  x*  ■  x  in  steps  of  Ax'  ■  Ax. 

In  both  cases,  simple  step  integration  is  used.  That  is 

/f(x')dx»  -  l  (Ax! ) [f (x{)  1  (64) 

all  i  x  1 
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( 
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The  integration  over  T'  is  carried  out  parallel  to  the 
integration  of  (56)  and  (57)  over  r  (using  space  as  a 
pseudo  retarded  time)  and  simultaneously  with  the  increase 
in  t  as  the  space  integrations  are  repeated  for  each  new  t. 

This  rather  involved  approach  to  solving  equation  (27) 
is  necessary  to  save  running  time.  A  direct  approach,  with 
separate  integrations,  would  at  least  triple  or  quadruple 
the  total  running  time  required  for  execution  of  the  code. 

Integration  of  the  Field  Equations 

For  each  t,  equations  (56)  and  (57)  are  integrated  from 
r  -  RMIN  tor*  RMAX  in  steps  of  Ar  *  (RMAX-RMIN) /NDELR 
using  the  4**1  order  Runge-Kutta  method.  Then  the 
magnitude  of  E  is  found  from  the  two  components  and  the 
result  is  stored  in  the  E  array,  t  is  increased  by  At 
and  the  whole  process  is  repeated  until  T  reaches  TMAX. 

On  completion  of  the  iterations,  each  member  of  the 
E  array  is  multiplied  by  RMAX/r  (equation  62).  Then 

w aTgcb 

the  E  array  is  searched  to  find  the  peak  value. 

Outputs 

There  are  several  output  options  available  in  the 
code.  The  basic  output  it: 

1.  Gamma  yield  and  altitude  of  burst. 

2.  Target  coordinates  from  ground  zero. 

3.  Distance  from  burst  to  target. 

4.  A  message  indicating  whether  the  direct  or  the 
reflected  wave  is  being  calculated. 
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5.  The  time  period  covered  by  the  calculation. 

6.  The  tine  when  the  peak  value  occured. 

7.  The  peak  value  of  E  at  the  target. 

8.  The  t  and  E  arrays. 

In  addition,  a  linear  and  a  log-log  plot  of  E(t)  can 
be  obtained.  Also,  a  listing  of  the  values  of  E  at  the 
bottom  of  the  absorption  region  for  each  t  can  be  obtained. 
Either  or  both  of  these  two  options  can  be  added  to  the 
basic  output. 


c 
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IV.  Results  and  Input  Parameter  Variation 
The  output  from  a  typical  run  is  shown  in  Pig.  4. 


The 

E(t)  calculated  during  the  run  is  shown  in  Pig. 

5. 

The 

input  data 

for  this  run  was: 

X 

■ 

0  meters 

(65a) 

Y 

■ 

0  meters 

(6Sb) 

Z 

■ 

0  meters 

(65c) 

HOB 

■ 

100  km 

(6Sd) 

\ 

U 

.001  kt 

(65e) 

Bo 

m 

2(10) ~5  wb/n2 

(6Sf) 

Dip  Angle 

m 

20* 

(6Sg) 

NDELR 

m 

SO 

(65h) 

TMAX 

m 

20  shakes 

(6Si) 

The  CDC  6600  Computer  required  191  sec  and  33000g  words  of 
central  memory  to  execute  this  run. 

The  peak  value  of  E(  6400  V/m,  obtained  in  this 
run  compares  favorably  with  Karzas-Latter's  order  of 
magnitude  estimate  of  104  V/m  (Ref  2)  from  similar  input 
data. 

In  order  to  gain  a  better  knowledge  of  the  operating 
capabilities  of  the  code,  the  effect  of  varying  input 
parameters  one  at  a  time  was  studied.  The  basic  set  of 
parameters  used  was: 

( ) 
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THE  BURST  WITH  GAMMA  YIELD  CF  1.000E-03  KILGTCNS 
IS  AT  AN  ALTITUCE  OF  1.000E4C2  KILOMETERS* 


Typical  Run 


EFIELD  (V/M)  «10* 

cp.OQ  20.00  40.00  60.00  80.00  10(1.00 

*•  i  i  i  i 


OO’CP 
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X 

■ 

0  meters 

(66a) 

Y 

■ 

0  meters 

(66b) 

Z 

■ 

0  meters 

(66c) 

HOB 

■ 

100  km 

(66d) 

Yv 

* 

.001  kt 

(66e) 

Each  of  the  above  parameters  was  systematically  varied 
while  holding  the  others  constant.  The  other  inputs  were 
held  constant  at  the  values  shown  in  equations  (65). 

The  results  of  the  variation  in  X  are  shown  in 
Pig.  6.  Since  the  X  axis  is  perpendicular  to  the  magnetic 
field  the  symmetry  about  X  ■  0  is  expected.  The  decrease 
in  peak  value  of  E  for  increasing  distance  from  ground  zero 
is  due  to  the  increasing  distance  from  the  burst. 

The  results  of  the  variation  in  Y  are  shown  in  Fig.  7. 
Here  the  peak  values  of  E  depend  on  the  angle  between 
r  and  1Qf  0.  When  0  •  180*  (A  -  -70*  and  Y  •  -275  km) 
the  peak  E  drops  to  zero.  The  maximum  peak  E  is  skewed 
toward  A  ■  20*  (0  ■  90*  and  Y  ■  36  km).  The  maximum  is 
not  exactly  at  A  ■  20*  because  of  the  increased  distance 
from  the  burst.  These  characteristics  are  expected  since 
an  electron  moving  perpendicular  to  the  magnetic  field  would 
feel  the  strongest  acceleration  frrm  it  while  an  electron 
moving  parallel  to  the  magnetic  field  would  feel  no  accel¬ 
eration  at  all. 

The  results  of  variation  in  Z  are  shown  in  Pig.  8. 

In  this  case,  both  the  direct  and  the  reflected  waves  were 
calculated  at  each  point  below  the  top  of  the  absorption 

32 
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region.  Note  that  Y  ■  100  km  for  these  runs.  As  expected, 
the  direct  wave  falls  off  rapidly  as  the  target  altitude 
passes  through  the  absorption  region,  since  less  of  the 
absorption  region  contributes  to  the  wave  with  each  increase 
in  altitude.  The  crossover  point  where  the  reflected  wave 
becomes  the  largest  occured  at  2S  km  in  this  case.  Above 
ground  zero  the  crossover  point  was  29.4  km.  The  altitude 
of  the  crossover  point  is  both  yield  and  geometry  dependent. 
It  is  necessary  for  the  user  to  calculate  both  waves  when¬ 
ever  there  is  any  doubt  which  one  is  the  largest. 

The  reflected  wave  calculation  assumes  100%  reflection 
from  the  ground  and  no  attenuation  in  the  absorption  region 
or  the  ionosphere.  These  assumptions  are  reasonable  if  it 
is  recalled  that  only  the  high  frequency  component  is  being 
considered  and  that  it  requires  at  least 

- —  «  133  p  sec  (67) 

3(10)°m/sec 

for  the  wave  to  leave  the  absorption  region,  reach  the 
earth,  be  reflected,  and  return  to  the  absorption  region. 

This  length  of  time  is  enough  for  a  significant  number  of 
the  free  electrons  to  recombine  and  reduce  the  effective 
conductivity  of  the  absorption  region. 

The  results  of  variation  in  HOB  are  shown  in  Fig.  9. 

For  all  values  of  HOB  attempted  below  60  km  the  code  went 
unstable.  Infinite  values  for  E  were  obtained  which  resulted 
in  abnormal  termination  of  the  calculations.  This  is 
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expected  since  the  burst  is  assumed  to  be  distant  from 
the  absorption  region  (equations  9  and  10). 

The  results  of  variation  in  gamma  yield  are  shown  in 
Pig.  10.  For  all  gamma  yields  attempted  above  60  tons 
the  code  went  unstable,  giving  infinite  values  for  E. 
However,  the  instability  always  occured  at  times  later 
than  the  natural  peak  value  of  E.  For  example,  with  80 
tons  of  gamma  yield,  the  natural  peak  occured  at  1  shake 
and  the  instability  occured  at  10  shakes.  By  using  the 
natural  peak  value  and  ignoring  the  instability,  reasonable 
values  for  peak  E  were  obtained  up  to  1  kt  of  gamma  yield. 


Peak  E  (Volts/meter) 
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V.  Discussion  and  Recommendations 


Limitations 

Most  of  the  limitations  of  the  code  are  inherent  in 
the  model  upon  which  it  is  based.  Approximations  such 
as  a  flat  earth,  a  uniform  magnetic  field,  and  constant 
speed  Compton  electrons  can  be  improved  only  by  changing 
the  model.  In  addition,  the  effect  of  the  self  generated 
electromagnetic  fields  on  the  motion  of  the  Compton 
electrons  is  ignored,  as  is  recombination  of  both  primary 
and  secondary  electrons.  The  possibility  of  a  single  gamma 
ray  interacting  to  produce  more  than  one  Compton  electron 
is  not  allowed.  In  the  absorption  region  the  contribution 
of  the  non-propagating  radial  component  of  the  electric 
field  is  neglected.  Also,  the  model  is  not  easily  adapted 
to  multi-group  gamma  transport,  or  to  multiple  burst 
calculations. 

The  code  calculates  only  the  effect  of  the  gamma 
rays.  The  user  must  keep  in  mind  that  X-ray  generated  EMP 
becomes  important  for  bursts  above  100  km. 

The  code  does  not  account  for  the  increase  in  altitude 
of  the  absorption  region  for  slant  angles  (angle  A  in 
Fig.  1)  greater  than  60*  which  is  indicated  by  Latter 
and  LeLevier  (Ref  4). 

Since  97%  of  the  running  time  of  the  code  is  used  for 
numerical  iteration  it  is  not  practical  to  adapt  the  code 
to  run  more  than  one  target  at  a  time.  Two  targets  would 
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merely  double  the  running  time,  so  it  is  simpler  to  just 
make  two  runs.  Typical  requirements  are  200  seconds 
running  time  with  33000g  words  of  central  memory  on  the 
CDC  6600  computer  using  NDELR  ■  50  and  TMAX  ■  20  shakes. 

Uses 

The  code  can  be  used  to  calculate  the  peak  value  of 
the  E  field  at  a  target,  anywhere  on  or  above  ground  level, 
resulting  from  a  nuclear  burst  above  60  km  altitude  with 
a  gamma  yield  up  to  60  tons.  Either  the  direct  or  the 
ground  reflected  wave  can  be  calculated.  With  special 
care,  bursts  up  to  1  kt  of  gamma  yield  can  be  used. 

Recommendations 

In  the  interest  of  accuracy,  the  targets  should  be 
located  such  that  the  slant  angle,  A,  is  between  -60*  and 
♦60*. 

By  accepting  a  much  longer  running  time  the  accuracy 
and  hopefully,  the  stability  of  the  code  could  be  improved 
by  using  a  smaller  step  size  in  the  integration  of  the 
Compton  current  equations.  Reducing  the  step  size  from 
one  tenth  of  the  Compton  lifetime  to  one  shake  would 
require  approximately  ten  times  as  much  running  time  as 
the  code  presently  requires.  This  possibility  should  be 
investigated  further  to  determine  the  optimum  step  size 
for  obtaining  the  best  relationship  between  accuracy  and 
running  time. 
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Another  possibility  for  increasing  the  accuracy  and 
stability  of  the  code  is  to  reduce  the  step  size  in  r.  The 
present  code  has  the  capability  of  dividing  the  absorption 
region  into  S00  steps  in  r  along  the  line  of  sight.  Of 
course,  the  running  time  required  for  500  steps  is  ten 
times  that  required  for  SO  steps.  A  modification  of  the 
code  to  allow  more  than  500  steps  would  increase  the  amount 
of  computer  core  required  as  well  as  increasing  the  running 
time.  This  provides  another  area  for  investigation  to 
determine  the  best  trade  off  point  between  accuracy  and 
running  cost. 

These  two  possibilities  could  be  investigated  with 
minor  modifications  to  the  present  code.  However,  the 
^  computer  time  required  would  be  considerable. 

In  addition,  there  are  numerous  possibilities  for 
improvements  in  the  model  itself.  Some  of  the  more 
important  ones  are; 

Using  multigroup  gamma  transport. 

Using  multigroup  Compton  electrons. 

Allowing  angular  distribution  of  Compton  electrons. 

Using  self  consistent  electromagnetic  fields. 

Including  the  low  frequency  components. 

Each  of  these  would  require  major  modifications  to  the 
present  code. 

P 
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EMP  Code  User's  Guide 
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EMP  Code  User's  Guide 

The  code  is  run  the  same  as  any  other  Fortran  Extended 
program,  but  due  to  the  running  time  it  should  be 
converted  to  binary  form  before  execution.  The  plotting 
subroutine  requires  an  on-line  plotter  and  both  linear 
and  log  plotting  library  subroutines. 

The  input  data  is  read  in  the  following  order: 

Data  card  #1,  using  FORMAT  (7F10.0,  215),  contains; 
X,Y,Z  The  target  coordinates  in  meters 
HOB  The  height  of  the  burst  in  kilometers 
(60  km  <  HOB) 

GAMYLD  The  gamma  yield  in  kilotons 
(GAMYLD  <  1  kt) 

BFIELD  The  Earth's  magnetic  field  in  wb/m2 
BANGLE  The  magnetic  field  dip  angle  in  degrees 
NDELR  The  number  of  steps  in  r  taken  through 

the  absorption  region  (50  <_  NDELR  <_  500) 
OUT  The  output  control  parameter 
Data  card  #2,  using  FORMAT  (13),  contains; 

ITER  The  time  period  covered  by  the  iterations 
in  shakes  (10  <  ITER  <  100)  (ITER  -  TMAX) 
Data  card  #3,  using  FORMAT  (4F10.0) ,  contains; 

A  Pomranning  constant  a  in  inverse  shakes 

B  Pomranning  constant  $  in  inverse  shakes 

RN  Pomranning  constant  N  in  shakes 

TO  Pomranning  constant  Tq  in  shakes 
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Default  values  are  provided  for  BANGLE*  BFIELD,  and 

2 

NDELR.  They  are  40**  0.00002  wb/ra  ,  and  50  respectively. 

If  these  default  values  are  desired*  zero  must  be  punched 
in  their  respective  card  fields. 

The  ground  reflected  wave  at  the  target  is  obtained 
by  reading  in  the  target  altitude*  Z *  as  a  negative  number. 
For  any  target  within  the  absorption  region*  both  the  direct 
and  the  ground  reflected  wave  should  be  calculated  to 
determine  which  one  is  the  strongest. 

For  values  of  GAMYLD  between  0.06  kilotons  and  1.0 
kilotons  the  code  will  most  likely  go  unstable.  This 
instability  occurs  after  the  real  peak  has  been  calculated* 
but  the  peak  value  printed  out  may  not  be  the  real  peak. 
Since  execution  is  terminated  when  the  field  becomes 
greater  than  1E15  V/m*  the  array  search  can  result  in  a 
false  peak  value.  In  this  case,  the  arry  itself  (or  the 
plot)  can  be  used  to  determine  the  real  peak  value. 

Increasing  NDELR  makes  the  step  size  in  r  through 
the  absorption  region  smaller  and  the  calculation  becomes 
more  accurate.  However*  total  running  time  varies  directly 
with  changes  in  NDELR.  For  example*  using  NDELR  •  100 
instead  of  NDELR  ■  50  will  approximately  double  the 
running  time  required  for  NDELR  ■  50. 

There  are  four  output  options  provided.  Option  0 
prints  out  the  informative  messages,  the  calculated  peak 
value  at  the  target,  the  E  array*  and  the  t  array.  Option 
1  adds  a  linear  plot  of  the  first  20  shakes  and  a  log-log 
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plot  of  100  shakes  of  E  as  a  function  of  t  at  the  target. 
Option  2  includes  both  Option  0  and  Option  1  and  adds  a 
printout  of  E  and  a  as  a  function  of  t  at  thd  bottom  of 
the  absorption  region.  Option  3  deletes  the  plots  from 
Option  2.  The  last  two  options  are  primarily  for  debugging 
since  a  partial  printout  is  made  for  each  completed 
iteration  even  if  execution  is  terminated  before  the 
iterations  are  completed.  The  first  two  options  are 
best  for  production  runs. 

The  only  requirements  on  the  Pomranning  constants 
are  N  must  be  chosen  such  that  equations  (61)  and  (62) 
are  satisfied,  all  of  then  must  be  positive,  and  o  >  6. 

Increasing  ITER  also  increases  the  running  time. 

For  ITER  •  10  shakes,  running  time  is  approximately  180 
seconds  on  the  CDC  6600  computer.  For  ITER  ■  100  shakes, 
running  time  is  approximately  340  seconds.  A  good 
compromise,  which  gives  nice  looking  plots,  is  ITER  »  20 
shakes  with  a  running  time  of  approximately  200  seconds. 

In  binary  form,  the  code  requires  33000g  words  of 
core  on  the  CDC  6600  computer. 
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EMP  Code  Flow  Charts 
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PROGRAM  CONTRL 


[START] 


^  IS  TARGE! 
ABOVE  49  kn 


?  IS  TARGEl 
JELOW  GROUND 


PRINT 

REFLECTED 

WAVE 

MESSAGE 


PRINT  DIRECT  WAVE  MESSAGE 
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I  CALCULATE  ANGLES  A  AND  0  1 

1  CALCULATE  RMIN  AND  RMAXl 

ICALL  EFIELD 

TO  GET  E(rJI 

ICALCULATE  FIELD  AT  TARGETl 

i - ' 

I  SEARCH  ARRAY  TO 

FIND  PEAK  VALUEl 

PRINT  THE  TIME 
PEAK  VALUE  OF  F 
TIME  AND  FI 

- 1 

PEAK  OCCURRED, 

'I  ELD  AND  THE 

ELD  ARRAYS 

?  IS  PLOT  DESIRED 


YES 


CALL 

ELGPLT 
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SUBROUTINE  EFIELD 


[START] 


I  READ 


IN1TALIZE  ARRAYS  AND  CONSTANTS 


INCREASE  T 


INITIALIZE 

INTERMEDIATE 

CONSTANTS 


CALL  COMPTN 

TO  GET  Jt,  J* 

0  (p 

ICALL  CONDC 

I  TO  GET  <x| 

• 

CALL  RNGKUT  TO  GET  NEW  E, 


CALL  RNGKUT  TO  GET  NEW  E_ 
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SUBROUTINE  COMPTN 


[START 


INITIALIZE  CONSTANTS 


CALCULATE  NEW  AND  NEW 

0  QP 


CALCULATE  NEW  PRIMARIES  FOR  NEGATIVE  T 


CALCULATE  NEW  PRIMARIES  FOR  POSITIVE  T 


INCREASE  T* 


IS  COMPTON  LIFETIME  ST*  ? 


CALCULATE  RATE  OF 
PRODUCTION  OF  SECONDARIES 


[RETURN' 
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c 

SUBROUTINE  CONDCT 
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SUBROUTINE  ELGPLT 


(ST 

\RT) 

PLOT  E(t)  ON  1 

LINEAR  SCALES  | 

CLIP  SMALL  VALUES  OF  E(r)  1 

PRINT  NEW  E(f)  ARRAY  | 

PLOT  E(t)  ON  LOG -LOG  SCALES  | 

• 

RETURN 
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SUBROUTINE  RNGKUT 


(START) 


UNITIALIZE  CONSTANTS  1 

I DEFINE  FUNCTION  1 

CALCULATE  NEW  E 
USING  RUNGE-KUTTA 

(RETURN) 
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EMP  Code  Listlni 
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